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Abstract. The problem of accurately determining the equation of state of nuclear and neutron 
matter at density near and beyond saturation is still an open challenge. In this paper we will 
review the most recent progress made by means of Quantum Monte Carlo calculations, which 
are at present the only ab-inito method capable to treat a sufficiently large number of particles 
to give meaningful estimates depending only on the choice of the nucleon-nucleon interaction. 
In particular, we will discuss the introduction of density-dependent interactions, the study of 
the temperature dependence of the equation of state, and the possibility of accurately studying 
the effect of the onset of hyperons by developing an accurate hyperon-nucleon and hyperon- 
nucleon-nucleon interaction. 



1. Introduction 

In the last decade the connection between astrophysical observations and the accurate 
determination of the inter-nucleon forces has become stricter and stricter. This is partly due to 
the availability of qualitatively much improved data on radii and masses of neutron stars, which 
now permit, in principle, to exclude some of the models that have been introduced over time. 
For instance, the recent discovery of very massive neutron star (with estimated mass > 2Mq) jl]), 
and the general agreement on the fact that maximum masses exceed the value IAMq, seem to 
exclude the hypothesis that at density of order 2/9o, with /?o = 0.16 fm~^, a significant fraction of 
hyperons appears, that, according to existing calculations, would give a maximum mass largely 
below the currently accepted limit OS]. However, our present knowledge of both many-nucleon 



interaction and of the hyperon-nucleon interaction is too scarce to give to these results any 
prediction meaning. One of the most important shortcomings of our present knowledge lies 
in the fact that calculations are performed using approximate interaction in conjunction with 
approximate methods, and this fact essentially prevents to draw definite conclusions. 

In this paper we will show some of the recent progress that has been made towards an ab-initio 
treatment of the many-nucleon problem at high densities, by means of Auxiliary Field Diffusion 
Monte Carlo (AFDMC) and Fermi Hyper-Netted Chain (FHNC) calculations. In particular, we 
will focus on the recent rediscovery of the idea of modeling many-body interactions by means of 
density-dependent potential [HIS], the inclusion of finite temperature effects, and the preliminary 
work on the inclusion of three-body hyperon-nucleon-nucleon forces in the treatment of a mixed 
hyperon/nucleon matter. 



2. Auxiliary Field Diffusion Monte Carlo 

We briefly recall some notions about the AFDMC method. It is part of a wider class of 
algorithms implementing a stochastic procedure for projecting an initial trial state of an arbitrary 
Hamiltonian onto the state of lowest energy belonging to the same subspace in the Hilbert space 
{(j)n} of the eigenstates of the Hamiltonian itself: 

hm r) = hm V c„e-^(^-^o = c^<^^. (i) 

n 

The index a represents a specific symmetry for the state described by a set of quantum numbers, 
and/or the symmetry imposed by the Pauli principle. The constant Eq is a reference energy 
needed to control the normalization of the propagated states. The propagation is obtained 
by expanding the initial state on a finite number eigenstates of the position and spin/isospin 
operators Each initial point is then propagated by sampling a kernel, which, in general, 

is approximated by a Trotter-Suzuki break up of the original propagator for small imaginary 
time-steps: 

(x'|e-^"^|x)(x|^') = (x'|e-^"[^-^o]|2;')(x'|e-^"^|x)(x|^) +o(At). (2) 

The resulting expression for the propagation of a set of configurations of A particles can be 
summarized in integral form: 



. r ^,^2 



^(^', Ar) = (y^^^) / di2exp[-Ar(y(i?') - E^)] exp 



{R-R' 



^{R,0). (3) 



Therefore, the state is propagated by sampling the propagator and applying it to a set of 
points representative of the initial wave function. The points are iteratively displaced and 
weighted according to the various terms in the propagator, until a sufficiently long imaginary 
time is reached. In general the distribution of points is importance sampled, in the sense that 
it is multiplied by some approximate expression of the wave function of the state searched in 
order to avoid the large fiuctuations in the normalization coming from the possibly divergent 
behavior of the potential energy. The scheme works flawlessly only if the projected state is the 
absolute ground state of H, i.e., if it is a function symmetric under particle exchange. For a 
many Fermion system, as in the case of nucleons, the calculation is proved to be unstable, in 
that the variance/average ratio diverges exponentially. In this case it is customary to apply 
some constraint on the sampled configurations in order to maintain the overlap between the 
sampled distribution of points and the state onto which we want to project. These procedure 
are approximate, because they introduce a (usually small) bias on the estimates by effectively 
changing the Hamiltonian used in the propagation. 



An additional source of difficulty is introduced by the fact that, for many-nucleon systems, the 
potential is usually non-local. Local expressions of the two-body interaction, like the Argonne 
AV18 potential, are written as a sum over operatorial components: 

18 

V{n,) = Y,Vp{nj)Op. (4) 
p=i 

For the simple AV6 case, discussed in this paper, the operators considered are: 

Op = (1, Ti ■ Tj) (g) (1, cTj • aj, Sij), (5) 

where Sij is the tensor operator. From the point of view of the Diffusion Monte Carlo (DMC) 
algorithm, the use of such potentials imply that the factor exp[—AT{V(R') — Eq)] cannot be 
treated as a local "weight" for the propagated points, but rather behaves as an imaginary time 
propagator for the spin/isospin degrees of freedom. In the AFDMC algorithm the problem 
is circumvented by recasting the interaction in a spin/isospin independent part Vsi and a 
spin/isospin dependent one V^^. The second is in turn written as a bilinear form in the 
spin/isospin operators. For simplicity, let us drop the isospin components. The interaction 
becomes: ^ 

V = Vsi + - ^ ai^^Ai^^j^sf^j^s- (6) 

By diagonalizing the matrix A, and indicating as and A„ the corresponding eigenvectors and 
eigenvalues, we have: 

^ 3A 

y = Vs^ + ^J2^-0l (7) 

n=l 

where: 

0„ = Y.^il'a,,s. (8) 

j,s 

At this point, in the small At limit, the spin dependent propagator can be recast by using 
the Hubbard-Stratonovich transformation in order to reduce the operatorial dependence from 
quadratic to linear: 

3A 3A „ 2 



n=l 71=1 

The main advantage of this procedure is in avoiding the necessity of a sum over all the 
components of the wave function corresponding to different two-nucleon states. This number 
grows exponentially with the number of nucleons, and prevents the straightforward application of 
Green's Function Monte Carlo (GFMC) techniques for systems larger than A = 12 [6j. However, 
the price to pay is the introduction of a set of auxiliary degrees of freedom by means of which 
rotation of the spinorial components of the single particle functions are sampled. Treating a 
sufficiently large number of nucleons in a periodic box is a necessary condition for studying 
the properties of infinite matter without suffering too heavy finite-size effects. A study of the 
convergence of the results to the thermodynamic limit can be found in Refs. [71 [8]. 

The AFDMC algorithm actually in use exploits a number of technicalities that are not 
reported here, including an extension of the importance sampling concept, and a particular 
form of constraint to avoid the problem of the exponential growth of the variance. All such 
details have been largely discussed in the existing literature. For neutron matter and neutron 
drops the method can be extended to Hamiltonians including spin-orbit (AV8') and three-body 
interactions |9l HOl IH] . In the case of nuclear matter the method is still limited to the use of 
the 6 operators mentioned above 
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Figure 1. Equation of state of symmetric and neutron matter from AFDMC calculations with 
a density dependent interaction (see legend). The curve computed at /3-equilibrium with and 
without inclusion of muons is also reported. 



3. Density Dependent Interactions 

An accurate description of the many-nucleon interaction in nuclear matter should in principle 
rely on the forces that describe the binding energies of nuclei. However, at present, the simple 
transposition of the Hamiltonian as is from nuclei to the infinite systems seems not to give 
satisfactory results. This is obviously due to the difficulty of developing a true ab-initio theory 
of the many-nucleon interaction. In order to be operative, and provide to astrophysicists some 
significant estimates of the equation of state, it is possible to implement a two-body potential that 
effectively includes many-body effects through a density dependence of the coupling constants. 
Following Lagaris and Pandharipande, the AV6 potential, can be recast by introducing a density 
dependent intermediate part: 

Vij = + e~"'^''vi + Vr- (10) 

In order to reproduce the saturation density and energy, it is then necessary to introduce a 
purely phenomenological attractive contribution which includes a symmetry term: 
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The three parameters are fitted, by means of AFDMC calculations, towards the saturation 
energy, saturation density and compressibility in nuclear matter jlj. Their values are 71 = 0.10 
fm^, 72 = —750 MeV-fm^, and 73 = 13.9 fm^ respectively. With these values is then possible to 
evaluate the energy per nucleon at densities above and below saturation. The resulting equation 



of state of the Symmetric Nuclear Matter (SNM) can be fitted by the expression: 



= Eo + Kp - pof + cip - pofe'^^^-^^\ (12) 

with Eo = -16.0 MeV, po = 0.16 frn'^, b = 520.0 MeV-fm^, c = -1297.4MeV-fm9, and 
7 = — 2.213fm^ [1]. In Fig. 1 we report the corresponding curve, together with the results 
of the equation of state of pure neutron matter (PNM) obtained computing the energy per 
neutron with the same density dependent interaction (DDI). Assuming a quadratic behavior 
in the proton/neutron number difference for the symmetry energy, it is possible to write an 
equation of state for a generic proton fraction Xp as follows: 

E{p,xp) = Esnm{p) + Cs (^^y (1 - 2xpf, (13) 

with Cs = 31.3MeV, and 7^ = 0.64. Typical values for these parameters have been quoted as 
« 31 - 33 MeV and 7^ w 0.55 - 0.69 by [13j and as Cs = 31.6 MeV and 7^ w 0.69 - 1.05 
by [II]. It should be noticed that, usually, the symmetry energy is constrained over a range of 
densities typical of nuclei, whereas here it was fitted the parameters over a very wide density 



range. This means that the parametrization of Eq. (13) should be accurate up to very high 
densities. 

In this way it is possible, density by density, to compute the equation of state at beta- 
equilibrium with electrons and muons, by assuming that the chemical potential of the leptons 
is described by that of an ultra-relativistic Fermi gas. The resulting curves are also reported in 
Fig. 1. 



4. Finite Temperature Corrections 

In the final stages of a supernova explosion, the dynamics of the resulting proto neutron star is 
driven by the equation of state at temperatures that initially reach values of about 20 MeV. In 
astrophysical models the temperature dependence of the equation of state is usually assumed 
to be very simple. So far, realistic calculations have been performed in the Brueckner-Hartree- 
Fock frameworkp^. In order to have a completely ab-initio description it would be necessary 
to turn to Path Integral based Monte Carlo methods that allow to compute expectations over 
the quantum thermal matrix. However, this step implies a number of technical difficulties that 
have not been faced yet. 

An interesting intermediate approach consists of estimating the thermal corrections to the 
zero temperature EoS by means of of a temperature dependent FHNC calculation [TB]. Such 
corrections can then be applied to the equation of state computed by AFDMC to obtain a result 
as accurate as possible. Several tests showed that by including the temperature effects using 
FHNC the contributions of the most important elementary diagrams cancel, and the difference 
is only weakly affected by the lack of such diagrams. 

By extending the variational chain summation to finite temperatures one has to face the so- 
called orthogonality corrections |16^ I17j. since wave functions used are not mutually orthogonal. 
The orthogonality corrections are not unique and, in general, their calculation requires the 
evaluation of the off-diagonal matrix elements of the Hamiltonian. At present no accurate 
method exists to efficiently compute such off-diagonal matrix elements. Recently it was proved 
that the orthogonality corrections to the free energy vanish in the thermodynamic limit |18| . In 
the spirit of Landau's theory of Fermi liquids, where it is assumed that the low- lying eigenstates 
of an interacting system have one-to-one correspondence with those of the noninteracting gas, 
we assume that a good approximation for the eigenstates of the interacting fermion system is 



given by the correlated basis states [TB] : 



[n^m = y ' ^ ^ = , (14) 

'^l[n^m{'S{n^<J^^J)) '^.[^.(k)] 

where are the single particle states of the non-interacting system . The nj(k) = 0,1 are 
the occupation numbers for single particle states labeled by k. The pair correlation operator 
defining the correlated wave function J^ij is taken to be 

6 

p=i 

where p = 1 — 6 operators are defined as in Eq. (5). Since the operators Ofj do not commute, the 
product of correlation operators is symmetrized with the symmetrization operator S to make 
the wave function antisymmetric. 

The upper bound for the free energy F(p, T) can be obtained by using the Gibbs-Bogoliubov 
variational principle 

F{p, T) < Fvip, T) = Tr {pyH) - TSv{p, T) . (16) 
pv is any arbitrary density matrix satisfying: 

Tipv = 1, (17) 

and Sv{p, T) = Tr {py Inpy) is the entropy derived from the density matrix py at temperature 
T. The equality in ( |16[ ) holds when py is the exact density matrix of the system. In order to 
obtain py it is necessary to start from an ansatz for a correlated effective Hamiltonian Hy, such 
that: 

- TVexp(-/3i?v) ' ^'^^ 

where /3 = l/T is the inverse temperature ( ks = 1). 

In practice Hy is chosen to be a one-body operator such that correlated basis states ( 14 ) are 
eigenstates of it: 



Hy\^i [ni(k)]) 



^ni(k)ey(k,p,r) 



l**K(k)]). (19) 



The eigenvalues of this Hy can be varied by changing the single-particle spectrum ey(k), 
which can be interpreted as the quasi-particle spectrum, and the eigenfunctions by varying 
the correlation operator Tij. 

With this choice of Hy, the calculation of the entropy Sy is trivial. At temperature T the 
average occupation number of a single-particle state is given by: 

^(k,p,r) = rr^— , (20) 

exp [p (e(k, p, T) - p[p, T))\+l 

where the chemical potential p{p, T) is required to satisfy the normalization condition: 



A = ^n(k,p,r), 

k 



(21) 



where A is the total number of particles in the system, and the entropy is given by: 

Sv{p, T) = -Y, P, T) In (n(k, p, T)) + (1 - n(k, p, T)) In (1 - n(k, p, T)) 
k 



(22) 



Since the correlated basis states (14) are not mutually orthogonal, the last equation is only 

an approximation if the variational Hamiltonian Hy is defined by Eq. (19). It is exact if only 

orthonormalized correlated basis states are be used. 

The calculation of Ev{p,T)/A = Tt(pvH) is very similar to the variational calculation of 

-2 _ 

ij 

^^^X^ = ^ + E diagrams(^;, i{r, p, T)) , 



the ground state expectation value Eq by expanding it in power of Ff^ — 1. Schematically, we 
have: 



2m 



(23) 



where m is the average bare mass of a nucleon and k^^ is a mean square momentum per particle, 



l^k2n(k,p,r). 



(24) 



The diagrams are many-body integrals involving the potential u, the correlation operator T and 
the finite temperature Slater function 



{r, ^) = T E exp(^kr)n(k, p, T) . 



(25) 



The single particle spectrum is parametrized as: 



21,2 



e(k,p,r) 



K^k 



2m 



l + A{p,T)exp{-B{p,T)k^) 



(26) 



It should be noted that in general it is possible for e(k, p, T) to have higher order terms in 
k. However, the free energy was found to be insensitive to any such dependence at moderate 
density range. Under this assumption, if i?(p, T) = 0, the spectrum is fully determined by an 
effective mass: 



^2 flde\ 



-1 



l + A{p,T), 



(27) 



m*(p, T) 

m m \k dk J 

that can vary with temperature and density. 

This scheme has been applied to compute the thermal corrections to SNM and PNM within 
the FHNC/SOC scheme. The EoS resulting by adding such corrections to the AFDMC result 
obtained with a DDI is plotted in Fig. 2. The following functional form provides a good 
parametrization of the numerical results in the required ranges of density and temperature: 



F(p, T, x)/A = Eip, x)/A + AFoip, T)/A + (1 - 2xfAFs{p, T)/A 



(28) 



-a 



X 



1/3 



+ (1 



a/3 



where a and /3 are almost independent on the isospin x. The fit is inspired by the Sommerfeld 
expansion, and resembles the excitation energy of a hot non-interacting Fermi gas |21] : 



{F - Fo) I A 



+ 0(T4 



(29) 
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Figure 2. Equation of state of Pure Neutron Matter (bottom panel) and Symmetric Nuclear 
Matter (top panel), computed by adding the estimate of the FHNC free energy estimation to 
the zero temperature results obtained by means of AFDMC calculations. Results are reported 
for temperatures ranging from to 30 MeV. The low density PNM EoS is fitted to the virial 
equation of state of Ref. [HI [20] 



For pure neutron matter at the normal density the parameter Oe has the value Sn"^ /{Sfip) 
0.03315 MeV-^ 



Other functions, entering the definition of (29), are the following: 
AFo{p,T)/A = 



AFsip,T)/A 



ax log p + 



m 


T + 







ct log^ p + (It ' 



(30) 
(31) 



The parameters have been fitted against the AFDMC results. Their values are the following: 




Figure 3. Entropy per nucleon in SNM (upper panel) and PNM (lower panel) as a function of 
the density of tlie system. 



OT = -0.15(2), bx = -0.38(4), ct = -0.008(1), dr = 0.06(3), er = -0.016(13), a = 0.047(23), 
(3 = 0.72(14) witli x^/n.d.f = 0.54. 

The entropy per nucleon S{p,T,x), which is a measure of thermal disorder, is calculated 



from the quasi-particle occupation probabilities n(k, /), T) using Eq. (22). It is also possible to 
compute S{p,T,x) by means of the following expression: 

5(/>,r,x) = -(^)^. (32) 

The values of the entropy computed by the two different procedures are in excellent agreement. 
This is a strong test of the quality of the variational calculation of F{p, T, x) as discussed in 
|22j . Notice that entropy production in multi- fragmentation events in heavy-ion collisions is a 
crucial quantity in the determining the mass fragment distribution. The entropy per particle is 
shown in Fig. [3] as a function of density and for various temperatures. The entropy increases 
with temperature, as physically reasonable, and decreases substantially with density. At low 
T, it is expected to approach a linear dependence due to the fact that, for a Fermi liquid, the 



relation between S and T should be approximately 



S»-iV(T = 0)T=^r, (33) 
in terms of the density of states at the Fermi surface. The specific heat is defined as 

Cy(p,r,x) = r(||)^. (34) 

The energy density e{p,T,x), sound velocity Cs{p,T,x) (in units of c) and adiabatic index 
T{p, T, x) are given by: 

e{p,T,x) = p [F{p,T,x)/A + {1 - x)mnc'^ + xnipC^) , (35) 

Cs{p,T,x) = (36) 

T{p,T,x) = ^cl (37) 

If r is constant, then the EOS becomes of the usual polytrope form, P ~ e'". 



5. Hyperons and hypernuclei 

The onset of degrees of freedom with strangeness S 7^ in nuclei and nuclear matter is a problem 
of great interest in the study of the properties of dense stars. At densities of about 2pQ, where 
po = 0.16 fm~^ is the saturation density of the nuclear medium, the chemical potential of the 
ultra- relativistic electron gas, determined by the /3-equilibrium condition becomes comparable 
with that of the T,~ hyperon, which becomes stable due to its larger mass and the consequent 
decrease in kinetic energy. This fact has strong consequences on the EoS, which becomes 
softer than that predicted by models in which strange degrees of freedom are absent. Present 
calculations of the equation of state (EOS) of dense matter including hyperons show that this 
softening leads to an unphysical limitation of the maximum observable mass of a star [2,j 3] . This 
seems to be an indication against the indirect evidence of the occurrence of this mechanism. On 
the other hand, most calculations neglect important pieces of the interaction, and in particular 
the hyperon-nucleon-nucleon (YNN) contribution, which, if repulsive on average, might in 
principle completely change this picture. Due to the fact that the YN interaction must be 
mediated by at least two pions, the YNN force appears at the same order, and cannot be 
assumed to be small. 

Recently we started an ambitious project, which should lead to a more accurate determination 
of the YN and YNN interactions combining AFDMC calculations and possibly new available 
data on separation energies in hyper nuclei. Due to the very limited availability of data on 
S-hypernuclei, it is necessary to focus on the AN and ANN interactions only. The starting 
model is that introduced by Usmani [23]. The system under study is a hypernucleus composed 
by A nucleons interacting through a two body AV6' potential, and one hyperon. We write the 
Hamiltonian of the system as: 

1=1 j>i=l i=l 

The explicit form of the A-nucleon potential is even in Ref. [23]. It essentially accounts for 
a two-pion exchange interaction. In principle, it should also include a contribution from Kaon 



exchange. This term contributes to the tensor components, but it also includes a space exchange 
term between the nucleon and hyperon degrees of freedom. As assumed, and partially justified, in 
other works like [24|l25j. we assume that the exchange term is quite negligible. The contribution 
can be considered as effectively included in the hyperon separation energy by the fitted value of 
the strength of the two-pion exchange term. 

The three-body YNN interaction is instead of the general form: 

Vynn = + V^^, (39) 

where is a dispersive term, an V"^"^ is once more a standard two-pion exchange contribution. 
It does not present any special difficulty in the AFDMC scheme, due to the assumed 
distinguishability of the A with respect to the nucleons. Particular care has to be taken in 
treating the center of mass contributions. 

The separation energy of the A particle is defined starting from the difference between the 
energies of the nuclear systems with and without the A hyperon: 

- B{K) = {Hn+a)a+a - {Hn)a , (40) 



The hypernucleus wave function is built starting from single particle orbitals computed by HF 
with a Skyrme I force. A general expression of the mean-field part of the wave function is then: 

I^A+a) = (jlhirM)^ ^f:i^^{rA)\^A) , (41) 

where 4>nijmi^^) orbital describing the hyperon. The function fhifKi) is a two-body scalar 

(Jastrow) correlation between the hyperon and a single nucleon and \^ a) is the correlated wave 
function describing the remaining A nucleons. In our calculation this function is defined as: 

1*^)= [\{ f^inj)^ <I>^,z(l,...,A), (42) 

where ^n,z is the Slater determinant of a set of single particle wave functions of neutrons 
and Z protons. Obviously, A = Z + N . 

The AFDMC calculations have been performed for different values of the parameters 
characterizing the hyperon-nucleon interaction. Despite the estimates of the A separation 
energies are still rather noisy, it was possible to determine a set of values of the potential 
parameters by which the experimental data are roughly reproduced at the same time for ^He 
and jv'O. With this set of parameters we gave a preliminary estimate of the separation energies 
by both including and excluding the YNN term in the Hamiltonian. Notice that we are making 
the strong assumption that the separation energy, being the difference of two terms, is not 
strongly influenced by the quality of the nucleon-nucleon Hamiltonian that we employ in our 
calculations. 

In Table 1 we report the results obtained in our simulations. As it can be seen, the first 
indication is that using a YN Hamiltonian only, the separation energy tends to be overestimated. 
On the other hand, the YNN term seems to give an overall repulsive contribution. This result 
would confirm the fact that the inclusion of three-body forces is necessary to correctly discuss 
the onset of hyperons in infinite matter, and their contribution to correcting quantities like the 
mass-radius relation, which are presently the subject of an open discussion. 



Table 1. Hyperon separation energy B{A) as computed by AFDMC simulations for a set of A- 
hypernuclei. Results are reported for Hamiltonians with a two body YN interaction, as described 
in the text, and for a Hamiltonian including a three-body YNN interaction. Experimental 
estimates are from Ref. |26| 
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